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Abstract 


This is the first in a two-part series of manuscripts describing numerical experiments on 
the influence of 2-30 km strip-like heterogeneity on wet and dry boundary layers coupled to 
the land surface. The strip-like heterogeneity is shown to dramatically alter the structure of 
the free-convective boundary layer by inducing significant organized circulations that modify 
turbulent statistics. The coupling with the land-surface modifies the circulations compared 
to previous studies using fixed surface forcing. Total boundary layer turbulence kinetic 
energy increases significantly for surface heterogeneity at scales between \jzi = 4 and 9, 
however entrainment rates for all cases are largely unaffected by the strip-like heterogeneity. 
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1. Introduction 

Land-atmosphere coupling is widely recognized as a crucial component of regional, conti- 
nental and global scale numerical models. Predictions from these large-scale models are 
sensitive to small-scale surface layer processes like heat and moisture fluxes at the air-soil- 
vegetation interface as well as boundary layer treatments (e.g., Garratt, 1993). The soil 
moisture boundary condition has a considerable influence on medium- to long-range weather 
forecasts and on simulated monthly mean climatic states {e.g., Rowntree and Bolton, 1983). 
Heterogeneous soil moisture conditions can occur at many scales both naturally (Mahrt et al., 
2001a) and through human modiflcation (Weaver and Avissar, 2001). Naturally occurring 
heterogeneity can arise through a number of mechanisms such as spatially varying soil or 
vegetation type and human-induced heterogeneity can result from spatially-varying agricul- 
tural practices that vary from farmer to farmer. Both types of heterogeneity can introduce 
dramatic variability in boundary layer surface forcing. 

There have been a number of prior investigations into the effects of non-homogeneous 
surface forcing on the PBL using high-resolution large eddy simulation (LES) {e.g., Hechtel 
et al., 1990; Hadfleld et al., 1991; Shen and Leclerc, 1995; Avissar and Schmidt, 1998; Al- 
bertson and Parlange, 1999; Albertson et al., 2001; Raasch and Harbusch, 2001; Esau and 
Lyons, 2002; Kustas and Albertson, 2003). Hechtel et al. (1990) used measurements to drive 
a relatively coarse resolution LES and found no significant modifications to homogeneous 
cases which they attributed to having a mean wind in their simulation. Hadfleld et al. 
(1991) studied relatively small-scale (1.5 and 4.5 km wavelength) strip-like heterogeneity 
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under free convective conditions and found that the heterogeneity induced organized mo- 
tions in the PBL. These motions increased horizontal and vertical velocity variances near 
the upper/lower boundaries and in the mid-PBL, respectively. Shen and Leclerc (1995) 
and Raasch and Harbusch (2001) showed that for small-scale patch-like heterogeneity, the 
intensity of the organized motions increased when the horizontal-scale of the heterogeneity 
increased to about the PBL height. When investigating land-surface forcing that is fixed in 
time, Avissax and Schmidt (1998) found larger-scale heterogeneity to continually intensify 
the organized motions when the heterogeneity increased from two to forty kilometers. Rather 
than specify the surface fluxes, which completely decouples the forcing from the overlying 
atmosphere, Albertson and Parlange (1999), Albertson et al. (2001), and Kustas and Al- 
bertson (2003) investigated land-surface heterogeneity using specified surface temperature 
and moisture. However, none of these studies used a dynamically coupled system where the 
soil responds to the atmosphere and vice versa. In addition, these studies are typically of 
relatively coarse resolution (12.5 m < Ax < 100 m, and 20 m < A 2 < 50 m) and, with 
the exception of Avissar and Schmidt (1998), focused on relatively small-scale heterogeneity 
(e.g., 250m < A < 5km). In our view, the shortcomings and sensitivities exhibited by large- 
scale numerical models are a consequence of both the inadequate modeling of the PBL and 
poor representation of the interaction with the land surface. In order to improve existing 
parameterizations, a more complete understanding of the mechanics and thermodynamics of 
air-soil interaction and the transport of water vapor by turbulent processes in the PBL are 
required. 
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This paper is the first of two manuscripts examining the interactions between the at- 
mosphere and the land surface using an LES model for the PEL coupled to a land surface 
model (LSM). Fine grids and large computational domains are used to examine the impact 
of a range of soil heterogeneity scales (A = 2-30 km) on PEL turbulence and the vertical 
flux of water vapor mixing ratio. The coupling between the PEL and the land-surface is 
found to be of fundamental importance in determining the PEL response. The heterogene- 
ity is shown to influence the PEL differently depending on the moisture content of the PEL 
air. In the second manuscript (Patton et al., 2003b), the organized motions induced by the 
heterogeneity are described in detail through phase-averaged statistics. 

2. Land- Atmosphere Coupling 

Over land surfaces during daytime hours, the surface energy balance (SEE) is thermodynam- 
ically coupled to the growth of the boundary layer through interactions between the surface 
sensible and latent heat fluxes, the atmospheric temperature and humidity and the growth 
rate of the boundary layer by entrainment (Raupach, 2000). Through the SEE, incoming 
solar radiation absorbed at the ground is partitioned to determine the energy available to 
the coupled land-atmosphere system. Key entities involved in the SEE, such as heat and 
water vapor fluxes, are influenced by PEL concentrations of these quantities which them- 
selves are influenced by their smface fluxes, forming a negative feedbaek loop (McNaughton 
and Raupach, 1996). 

The main focus of this research is on the PEL, but in an attempt to incorporate these 
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important feedbacks in the study, we coupled our well-established 3D, time-dependent, LES 
code (Sullivan et ah, 1996) to the NOAH (National Center for Environmental Prediction / 
Oregon State University / Air Force / Office of Hydrology) LSM, version 2.0 (Chang et al., 
1999). Numerous models of varying sophistication and empiricism have been proposed to 
describe the physics of soil- water- vegetation systems. Our choice of the NOAH model is 
based on the following desirable features: 1) the straightforward manner with which it can 
be coupled with a parent atmospheric model, 2) only a few empirical constants are needed 
to describe the land surface, 3) its good predictive capabilities, and 4) the robustness of the 
model itself. 

The LES predicts three-dimensional and time-dependent velocity fields {u, v, w) by nu- 
merically integrating a filtered set of incompressible Navier-Stokes equations under the 
Boussinesq approximation. To incorporate buoyancy effects, the LES also integrates a ther- 
modynamic energy equation (6) and an equation for water vapor mixing ratio (q), which 
are combined to calculate virtual potential temperature (9y = 9 -f O.610g). The effects of 
unresolved motions in the LES are estimated by following an equation for subgrid-scale 
(SGS) energy (e). In the soil, the LSM (NOAH) predicts vertical profiles of temperature 
(T) and moisture (0), by integrating a 1-D set of the soil thermodynamic and hydrologic 
equations (Mahrt and Pan, 1984; Pan and Mahrt, 1987). At every {x,y) grid point, the lower 
boundary conditions for the LES and the upper boundary conditions for the LSM define the 
coupling between the PBL and the land-surface. Atmospheric velocity boundary conditions 
are determined through dynamically calculating a drag coefficient from a prescribed surface 
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roughness length and the wind field at the first LES grid level, while surface temperature 
and moisture boundary conditions are determined dynamically at every grid point by solving 
the SEE. See Appendix A for other numerical details of the code. 

The SEE for an infinitesimally thin layer at the ground surface can be written (following 
Erutsaert, 1982; Chang et ah, 1999): 

Rn — H LE "i" G 4" Misc (1) 

where radiation fluxes toward the surface are taken as positive and all other fluxes are 
considered positive when leaving the surface. This sign convention is loosely based on the 
premise that radiation drives the energy of the surface and that the other energy fluxes are 
the means by which the surface disposes of the absorbed radiation. In Equation (1), Rn is 
the net radiation, which can be broken down into the sum of its parts: Rn — Rs{l — a) + 
^Rdn — where Rg is the incoming solar radiation, a is the albedo of the surface, e is 
the surface emissivity, Rdn is the downwelling long-wave radiation, a is the Stefan-Boltzman 
constant, and Tg is the surface temperature. H is the sensible heat flux, parameterized as 
—pCpChU{0a — Og) where p is the air density, Cp is the specific heat of the air, Ch is an 
exchange coefficient for heat, U is the wind speed immediately above the surface, and Oa and 
dg are the potential temperature of the air and ground surface respectively. LE is the latent 
heat flux, parameterized as —pLCeU{qa — Qs) where L is the latent heat of vaporization, 
Cg is an exchange coefficient for vapor, and and qg are water vapor mixing ratio of the 
air and ground surface respectively. G is the soil heat flux which can be parameterized 
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as kriQi) , where ^^(Oi) is the soil thermal conductivity that depends on the soil 
moisture 0i at the first soil level Z\ and soil type. Misc. represents other terms such 
as contributions from precipitation, melting snow, horizontal advection, or storage, which 
we ignore. Dynamically, the atmosphere and land-surface are coupled through the surface 
momentum flux. The surface momentum flux is parameterized as r = pul — pCmU^, where 
M* is the friction velocity and Cm is the drag coefficient. Through dynamic calculation of r, 
Rn, H, LE, and G the atmosphere is therefore directly coupled to the land-surface. 

This coupling is also a natural mechanism to undertake an investigation into the influence 
of land-surface heterogeneity on PBL structure. It is known that under free convective (no 
wind) conditions, imposed land-surface heterogeneity induces organized atmospheric motions 
{e.g., Hadfleld et al., 1991), but how these organized motions feed back to dry/cool the soil, 
and how that drying/cooling influences the organized motions has yet to be established. 
In the horizontal, landscape heterogeneity has been shown to influence the PBL structure 
at scales up to 40 km (Avissar and Schmidt, 1998). However, the typical (x,y,z) domain 
of LES studies of the PBL to date is about 5x5x2 km^ (e.g., Moeng, 1984; Hadfleld 
et al., 1991). Therefore, in order to incorporate heterogeneous land-surface forcing in our 
simulations, it is imperative that we consider computational domains and time-scales that 
are large compared to those routinely considered in LES of the PBL. Therefore, an LES 
that utilizes adequate resolution to capture turbulence and at the same time covers a large 
horizontal domain requires large numbers of grid points. For example, in a 30 x 5 x 3 km^ 
domain, 9 x 10® grid points are required for a mesh size of 50 x 50 x 20 m^. 
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3. The LSM: Deployment, Input Parameters, and Initial Condi- 
tions 

The LSM is implemented at every (x,y) grid point in the LES. Each of these grid points 
is covered by perennial grasses, with a roughness length for momentum Zo^ of 0.1 m, and 
overlays four silty clay loam soil layers at depths of [0.05, 0.20, 0.60, 1.00] m. Since the surface 
is covered by grasses, we assume a constant surface albedo of 0.2 across the domain for all 
cases. All snow and ice components of the LSM are turned off. The exchange coefficients for 
momentum, heat, and moisture {Cm, Ch, Cg) are evaluated locally at every timestep based 
on Paulson’s (1970) similarity functions. Ch and Cg are assumed equal. 

In the LSM, the parameter Cz defines the relationship between the roughness height for 
momentum Zo„ and the roughness height for heat Zo^: 

'' “* ( 2 ) 

To determine the value of Cg, following Brutsaert (1982) we assume that the exchange coeffi- 
cients for momentum (Cm) and heat (Ch) are related by the expression: CmjCh ~ 3, which 
is valid for bluff bodies (see Figure 5.4, Brutsaert, 1982). Assuming the drag law to be ap- 
plicable, we know that Cm — u'l/U'^ and Ch = 0^./{A6 U). Therefore CmjCh = f’’n{zo„/zo^), 
or Zo„/zo^ = e^. Substituting into Equation (2) gives Cz = 0.1. This is different from the 
default value (cz=0.2) in NOAH thereby increasing the efficiency of surface scalar exchange 
compared to NOAH’s default, but is derived directly from an understanding of land-surface 
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exchange at the scale of interest. 

In addition to modifying Czt the LSM needed another modification to be applicable for 
the relatively small-scale application discussed here. As it was obtained from the National 
Centers for Environmental Prediction (NCEP), the LSM assumes the convective scaling 
velocity [w, = zi (g^-|-0.61E))^/^] is the minimum wind speed to avoid singular solutions 
under zero wind (free convection) conditions. While this assumption may be reasonable for 
the larger-scale flows where the PBL is not resolved, «;* is a bulk PBL parameter and would 
therefore decouple the soil from the local overlying atmosphere. To remain consistent with 
our adopted coupling method (through drag law assumptions), we assume a local drag law 
formulation to determine a small convective velocity scale following Schmidt and Schumann 
(1989) such that the minimum wind speed at the ground surface is dynamically calculated 
from the following: 


«min = 0.07(|- A2CftA0„)5 (3) 

Uq 

where, g is the gravitational acceleration, Oo is a reference temperature, is the distance 
between the ground and the first grid point above the surface, Ch is the exchange coefficient 
for heat derived by the LSM and A0„ is the virtual potential temperature gradient between 
the same two heights. 


[Figure 1 about here.] 

Initial soil conditions were taken from ground-based measurements at the Little Washita 
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site on Day 193 of the Southern Great Plains 1997 (SGP97) experiment. This day and 
location were chosen because nearly an inch of rain fell on the preceding 1.5 days and they 
coincided with an intensive observing period. The observations were interpolated to the four 
soil levels in the LSM. Before using these as initial conditions for the LES-LSM runs, the 
soil model was run off-line for four complete diurnal cycles to allow the LSM to equilibrate. 
The noon-time conditions after these four diurnal cycles were used as the wet initial soil 
conditions (Figure 1). It is important to note that these wet initial soil conditions are such 
that the evaporation will be determined by a combination of the available soil moisture and 
the atmospheric demand rather than at the evaporative potential rate. Dry soil conditions 
were obtained by running the off-line LSM further for another seven diurnal cycles allowing 
the soil to dry under realistic forcing. Noon-time conditions after the seven day dry down 
period were chosen as dry conditions (Figure 1). Average initial soil conditions were obtained 
by picking the noon-time soil conditions that provided surface fluxes nearly equal to the 
noon-time fluxes averaged every day over the entire one-week dry-down period. 

4. Simulations 

The coupled LES-LSM simulations employ (600,100,144) grid points in the (x,y,z) directions 
representing a (30,5,2.88) km domain. This leads to a grid spacing of (50,50,20) meters in 
each of the (x,y,z) directions. 

A geostrophic wind of 0 m s~Ms imposed in both x and y directions, therefore the 
simulations are in the free convection limit. The only external forcing imposed on the 
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system is specified through the incoming solar radiation, which is set to 700 W m“^ and is 
constant in time for all cases. 

For all of the simulations, the initial potential temperature is constant with height (300K) 
below 790m. A capping inversion is imposed (3K/0.1km) between 790 m and 890 m, and 
the middle of this layer we call the initial inversion height, Z( = 840 m; above 890 m the 
stratification is 3 K/km (Figure 2). 

[Figure 2 about here.] 

Two different PBLs are simulated with varying moisture content. In the first set of cases, 
the PBL was initialized dry (1 g kg“^) throughout the entire domain (DP, dry PBL). The 
second set of cases was initialized with a relatively wet atmosphere (8 g kg“^) in the PBL, 
dropping sharply to 1 g kg~* across the initial inversion layer, between 790 and 840 m (WP, 
wet PBL). See Figure 2 for a visual representation of the initial DP and WP conditions. 

For each of the DP and WP cases, nine simulations were performed. Three of these are 
horizontally homogeneous cases with uniform soil conditions that are wet (WS), dry (DS), 
and average (AS) respectively. The subsequent six cases are initialized with horizontally 
heterogeneous soil conditions. The heterogeneity that is imposed occurs solely in the x- 
direction as a step-function change between wet and dry soil. We define A as the wavelength 
of one complete wet and dry cycle. The wavelengths vary as A = [2, 3, 5, 7.5, 15, 30] km. At 
initial times, the x-extent of a single patch (A/2) ranges therefore from 1.2zi to 17 . 82 ;j. The 
soil conditions are homogeneous in the ^/-direction for all simulations. 

The nomenclature used to discuss the cases from here on will adhere the following format: 
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The case name includes a description of the initial PBL moisture content followed by a 
description of the underlying soil conditions (separated by a hyphen). For example, case 
DP-WS refers to the run which was initiated with a dry PBL and a horizontally homogeneous 
wet soil condition, and case WP-7.5S refers to the case initialized with a wet PBL and a A 
= 7.5 km wavelength heterogeneous soil condition. The rest of the names follow accordingly 
(see Table 1). 

Each simulation is integrated out in time until the turbulence is in equilibrium with its 
external forcing. This time is defined by the necessary time for the boundary layer averaged 
turbulence kinetic energy and the boundary layer growth rate to become nearly constant. 
These criteria were typically satisfied within (2.8-3.3, 1. 9-2.5) hours for the (DP, WP) cases, 
respectively. The time steps for these simulations were variable but ranged between (0.7- 1.6, 
1.4-1. 7) s. Therefore, the simulations were integrated for about (6000, 4000) time steps before 
any averaging began. Statistics were calculated over the subsequent 4000 time steps, or about 
(4-9, 8-10) turnover times, where a turnover time is defined as For normalization 

purposes, we define 0, == Hf(pCpWt), 0* = E/{pw») and = (^ + O.610o E)/(pWt). See 
Table 1 for the horizontally- and time-averaged values of these quantities. 

[Table 1 about here.] 


5. Verification 

As mentioned above, large numbers of grid points are required to investigate desired scale 
range of heterogeneity while maintaining the ability to resolve the response of important 
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PBL scales of motion. Therefore, before performing the land-surface coupling, the NCAR 
LES code was ported from a shared memory (Cray vector supercomputer) environment to 
a distributed memory massively parallel {e.g., IBM SP3 supercomputer) environment using 
the Message Passing Interface (MPI). A single MPI decomposition in the vertical direction 
is used, which is natural since finite differences are employed in the vertical direction and 
a spectral representation in the horizontal directions. Initial tests of the new code indicate 
that the code scales well, i.e., the wall clock time decreases linearly as the number of MPI 
processes increases (not shown). Each simulation discussed in this manuscript took about 
(54,45) wall clock hours for the (DP, WP) cases on 48 IBM SP3 CPUs, or about (2592, 2160) 
GPU hours. In order to achieve acceptable efficiency in the coupled LES-LSM, it was also 
necessary to parallelize the LSM portion of the code across all MPI tasks. 

5.1 Comparison to Previous Studies 

Before embarking towards an extensive study with the newly developed coupled LES-LSM, 
tests were performed to verify that results generated were comparable to previous research. 
The investigation presented by Schmidt and Schumann (1989) is considered the first to exam- 
ine the free convective PBL in a robust and concise manner. The same case as Schmidt and 
Schumann (1989) was run with the new code. The main difference between the simulations 
is that the current simulation uses a coupled land-surface boundary condition. Through 1-D 
off-line LSM testing we devised a set of soil conditions and incoming solar radiation that 
would provide surface fluxes similar to those imposed in Schmidt and Schumann (1989). 
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Measurements from the Willis and Deardorff (1974) tank experiments provide observational 
support. 


[Figure 3 about here.] 

The variances and SGS energy from the coupled code are slightly larger in magnitude 
through the majority of the PBL than those of Schmidt and Schumann (1989) (Figure 3), 
which is consistent with the comparison of these codes presented in Nieuwstadt et al. (1993). 
Notable, however, is that the current coupled code outperforms Schmidt and Schumann 
(1989) code at reproducing the turbulent statistics from the Willis and DeardorflP (1974) 
tank experiments. These results provide confidence that incorporating the LSM at the lower 
boundary and porting the code to massively parallel machines has not dramatically altered 
the performance of the code. 

5.2 Evaluation of lateral domain size 

The imposed heterogeneity varies only in the z-direction, not in y. Therefore, we utilize 
a large domain in the avdirection and a relatively small domain in the ^direction. It is 
prudent to evaluate whether the relatively short lateral extent of the domain influences the 
flow that ensues. If the j^domain is too small, the organized motions that develop will 
become artificially coupled through the periodic boundary conditions. To evaluate whether 
our domain choice is large enough to represent that of an infinite domain, two additional 
simulations were performed. These cases are identical to case DP-15S except that the y- 
domain is either halved or doubled. Therefore, the two cases utilize 50 or 200 nodes to 
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represent 2.5 or 10 km in y, respectively. We chose to evaluate the DP-15S case since the 
heterogeneity strongly influenced the statistics (to be shown in Section 6.). 

[Figure 4 about here.] 

Following the analysis of Komminaho et al. (1996), too small a domain should reveal 
itself through amplification of two-point horizontal velocity correlations [i?uu(A^)] with in- 
creasing lag in the lateral direction (Ay). Comparison of R^u (Figure 4a) suggests that 
halving the domain size amplifies the lateral correlation of horizontal velocity (most notice- 
able at Ay = 1.2 km). Doubling the size of the lateral dimension has some impact on the 
correlations, slightly increasing the amplitude at Ay = 1.2 km and minimally decreasing 
the lateral distance corresponding to the first zero-crossing. Vertical profiles of normalized 
horizontal root-mean-square velocities (urms/^*) where Urms = {u"'^ + confirm that 

the extent of the lateral domain has little influence on predicted turbulence statistics (Figure 
4b). Therefore, through this analysis we are reasonably confident that our domain choice of 
30x5x2.88 km^ does not adversely influence the quantitative results and adequately repre- 
sents simulations of heterogeneity that infinitely span the lateral dimension. 

6. Results 

6.1 Instantaneous Fields 

6.1.1 Horizontally Varying Land-Surface Properties and Forcing 

[Figure 5 about here.] 


14 



I 


As an example, instantaneous horizontal slices of land surface properties from the last 
time step of simulation DP-AS (Figure 5) show evidence of the coupling between the PEL 
and the underlying soil. The surface heat fluxes (sensible, latent and soil) contain the 
signature of cellular convection suggesting that the exchange at the ground surface is driven 
by atmospheric demand. Also notice that the soil temperature and moisture at the first 
model level in the LSM provide evidence of spatially varied drying and warming; at this last 
time step, the soil has dried considerably since initialization as it began with a horizontally 
homogeneous distribution of 21.5% volumetric soil moisture. The soil temperature at the 
first model level has also increased from its initial horizontally homogeneous value of about 
309 K. Cases with heterogeneity reveal a similar evolution and pattern but those patterns are 
dominated by the disparity between wet and dry soil. Patton et al. (2003b) presents results 
and discussion of the surface fluxes as influenced by the heterogeneity-induced organized 
motions. 

6.1.2 The PBL 

[Figure 6 about here.] 

[Figure 7 about here.] 

For the DP cases, the influence of the varying soil moisture is evident in the flow visu- 
alization of vertical velocity and water vapor mixing ratio contours shown in Figures 6 and 
7. These figures show how the presence or absence of soil moisture variation dramatically 
alters the turbulent motions in the PBL. Over the region of high soil moisture content {e.g., 
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x= [0,3.75], [7.5,11.25], [15,18.75], [22.5,26.25] km. Figure 7) the convective plumes tend to 
be less vigorous compared to their counterparts in the region of low soil moisture content. 

At the same time, the shape of the PBL interface layer {i.e., the entrainment zone around 
2 : = 1.3 to 1.5 km) is also altered by soil moisture heterogeneity. Entrainment of low water 
vapor mixing ratio is observed over the region of high soil moisture and entrained dry air 
extends all the way to the ground (compare the second panels of Figures 6 and 7) In general, 
we see evidence of entrainment of dry air near the plume valleys, similar to that reported by 
Sullivan et al. (1998). The average normalized mixing ratio in the domain is about the same 
for the two cases (not shown), but in the DP-7.5S case, averaged vertically, the lowest PBL 
mixing ratio exists over the soil with the highest volumetric soil moisture; results presented 
in Patton et al. (2003b) will quantify the organized motions responsible for these patterns. 

The bottom two panels of Figures 6 and 7 show the r-variation of the associated instan- 
taneous surface fluxes and volumetric soil moisture at the uppermost LSM level. In contrast 
to studies with specified surface fluxes, the fluxes in the current simulations are spatially 
varying and respond to both the PBL and to soil temperature/moisture. The partitioning of 
energy by the LSM is evident in both figures, but especially in Figure 7 where over the wet 
soil a greater portion of the available energy goes into evaporating water and hence larger 
latent and smaller sensible heat fluxes compared to that over dry soil. Figures 6 and 7 are 
presented as examples; the picture is similar for the other heterogeneous cases. 
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6.2 Bulk Boundary Layer Statistics 
6.2.1 Surface Fluxes 

Since the only external input of energy to the simulations is the specified incoming solar 
radiation, the surface buoyancy forcing is dynamically determined through the LES- 
LSM coupling. Therefore, it is prudent to quantify the bulk forcing on the PBL. To do so, 
the horizontally- and time-averaged sensible heat flux {H), latent heat flux {LE) and the 
Bowen ratio {(3 = H/LE) are calculated. These quantities are presented -in Table 1. Note 
that and 0* can be calculated from and LE, and can be combined into the surface 
virtual temperature flux using = 6* + (0.61 ©*), where Oo is a reference potential 

temperature. is used throughout the manuscript as the scaling parameter for virtual 
potential temperature. Table 1 also includes the horizontally- and time-averaged boundary 
layer heights (zi) and the convective velocity (w«). Here, Zi is determined by finding the 
local height of the largest gradient in virtual potential temperature profile at every (x,y) 
location and averaging across all x and y. This method to determine Zi was first described 
as the gradient method in Sullivan et al. (1998), and further described as a wavelet method 
of minimum dilation in Davis et al. (2000). 

As expected, the horizontally homogeneous dry soil cases (DS) partition 7.S-7.7 times 
more incoming solar radiation into sensible heat compared to latent heat than does the wet 
soil case (WS). Also apparent from Table 1 is that our attempt to pick soil conditions for 
the average (AS) cases that resulted in average surface fluxes equal to the average of the 
WS and DS cases was not completely successful, e.g., of the AS cases is larger than the 
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average of those of WS and DS. Nonetheless, we can still use these AS cases as a point of 
reference. Scaling parameters Zi and reveal the expected influence of the magnitude of 
the surface buoyancy forcing 0^,*; the larger the larger are Zi and w^. Worth noting is 
that the land-surface model responds as expected to the moisture-state of the PBL; when 
the PBL is of higher average mixing ratio there is less atmospheric demand for moisture at 
the ground which is reflected in lower evaporative fluxes in the WP cases compared to the 
DP cases. 

6.2.2 Surface Exchange Coefficients 

[Figure 8 about here.] 

As mentioned earlier, at every {x,y) location in the domain the surface fluxes of momen- 
tum, heat and moisture in these simulations are internally determined through a coupling 
between a partitioning of the available energy energy via the SEB and a drag law. The 
efficiency of the transfer between the land-surface and the PBL is determined by drag coef- 
ficients {Cm, Ch) which are dynamically calculated from surface similarity theory assumed 
valid locally between the ground and the first LES grid level (which due to the staggered 
vertical grid system in the LES is 10 m for the velocities). Again, we assume that the 
exchange coefficient for moisture is the same as that for heat. 

On average, small-scale land-surface heterogeneity {\/zi ~ 1-2) increases the exchange 
coefficients for momentum and heat (Figure 8), while larger-scale heterogeneity {\/zi ~ 
4-9) tends to reduce these exchange coefficients. For flows encountering even larger-scale 
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heterogeneity {X/zi ~ 17), the coefficients decrease compared to \/zi ~ 4-9 heterogeneity 
and are tending back toward the values determined from the homogeneous cases. Cm and 
Ch for the X/zi ~ 4-9 cases are about 6% and 9% lower compared to the AS cases. It is 
understood that Cm and Ch are wind speed dependent (Schlichting, 1979; Brutsaert, 1982; 
Mahrt et al., 2001b), and that they decrease with increasing wind speed. The organized 
motions that develop as a result of the heterogeneity (see Patton et al., 2003b), are most 
intense for the A/zj ~ 4-9 cases. Therefore, the decrease in Cm and Ch for these cases likely 
results from a persistent increase in surface winds associated with the organized motions. 

Comparing the left and right panels of Figure 8, there is little or no influence of the initial 
PBL water vapor mixing ratio on the efficiency of the surface exchange which is consistent 
with Schlichting’s (1979) dimensional analysis suggesting that Cm should only be a function 
of Reynold’s numberh While the X/zi ~ 4-9 cases show the smallest Cm and Ch, their 
surface fluxes of heat and moisture are larger than the others (Table 1). This indicates that 
the gradients of heat and moisture between the ground surface and the first model level are 
larger on average for these cases than for the other heterogeneous soil cases. This increased 
gradient manifests itself through slightly lower average soil temperatures and volumetric 
moisture content combined with slightly higher PBL temperatures and water vapor mixing 
ratios (not shown). Although the effects are relatively small, this analysis implies greater 
transfer of heat and moisture from the soil to the PBL in the presence of X/zi ~ 4-9 scale 
heterogeneity. 

^Note that Ch depends on Cm (Brutsaert, 1982). 
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6.2.3 Turhvlence Kinetic Energy 


[Figure 9 about here.] 

Total turbulence kinetic energy (averaged over the PBL and time, ( )pBL,t) responds to 
land-surface heterogeneity (Figure 9). Heterogeneity at scales larger than 3 km increases 
the volume- and time-averaged turbulent kinetic energy. The increase in turbulence kinetic 
energy is due to the generation of organized motions (discussed in Patton et al., 2003b). 
Hadfield et al. (1991) also discussed circulations such as these in their LES simulations, but 
suggested that the circulation is weak since it only accounts for wl% of the kinetic energy 
in the PBL. The Hadfield et al. (1991) study, however, only looked at two wavelengths of 
heterogeneity, A = (1.5, 4.5) km, both of which fall at the small-scale end of this study and 
are in agreement with these results. As the scale of the heterogeneity increases towards a 
wavelength of 15 km, the circulation becomes stronger, and then decreases again at even 
larger scales. The TKE is most affected for cases DP-15S and WP-15S {X/zi ~ 9) which 
reveal an increase of (14, 19)% respectively compared to the horizontally homogeneous cases. 

6 . 2.4 PBL Height and Entrainment Rates 

[Figure 10 about here.] 

The growth of the convective layer is tied to the mixing process that occurs at the density 
interface between the well mixed PBL and the non-turbulent overlying stable layer (Sullivan 
et al., 1998). Apparent from Figures 6 and 7, under free convection, heterogeneous soil 
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moisture induces organized PBL motions that alter how the PBL operates. Therefore, it is 
anticipated that heterogeneity could modify the growth rate of the PBL. 

Figure 10 depicts the time evolution of the horizontally averaged boundary layer height 
for all cases. The times presented represent the entire averaging period for our statistics 
representing (4-9, 8-10) large eddy turnover times for the (DP, WP) cases, respec- 

tively. The average boundary layer growth is linear with time suggesting that the flow is 
quasi-steady. All heterogeneous cases fall between the horizontally homogeneous extreme 
cases (DS and WS). 

[Figure 11 about here.] 

Entrainment rates = dzi/dt) are calculated by performing a least-squares fit to the 
data presented in Figure 10, taking the slope of the fitted line to be We- Heterogeneity and 
the initial PBL moisture only slightly modify the normalized entrainment rates (Figure 11). 
Heterogeneity-scales on the order of the PBL depth tend to decrease normalized entrainment 
rates compared to the AS cases, although this decrease is only (5.6, 5.5)% for the (DP, WP) 
cases respectively. Heterogeneity-scales of (9, 5) times Zi slightly increase the entrainment 
rates for the (DP, WP) cases, but this increase is only (5.3, 3.8)% compared to the AS 
cases, respectively. These entrainment rates are also slightly higher than the homogeneous 
dry soil cases (DS). Even larger heterogeneity {X/zi ~ 18) tends to slightly decrease en- 
trainment rates. These results are counter to Avissar and Schmidt’s (1998) comments that 
heterogeneously-induced organized motions dramatically increase the rate at which the PBL 
entrains. However, they made these comments without actually calculating We/w^. Simi- 
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lar rates of entrainment are found for the homogeneous cases when comparing dry versus 
wet PBLs (compare the left and right panels of Figure 11). For both homogeneously- and 
heterogeneously-forced PBLs, the normalized entrainment rates presented in Figure 11 are 
noticeably higher (~ 0.04) than those presented in Sullivan et al. (1998) (~ 0.02) due to both 
the relatively weak inversion and to the strong surface forcing in these simulations compared 
to those in Sullivan et al. (1998). 

The combined results presented in Figures 9 and 11 indicates that the entrainment rate 
is not directly proportional to TKE, rather entrainment appears more dependent on the rate 
of heating. We can only make this statement as it applies to free convective PBLs interacting 
with and in equilibrium with strip-like heterogeneity. 

6.3 Horizontally Averaged Statistics 
6.3.1 Velocity Variances 

[Figure 12 about here.] 

[Figure 13 about here.] 

The moisture state of the PBL does not vastly influence horizontal and vertical velocity 
variances (Figures 12 and 13). Compared to the horizontally homogeneous cases, {u"^)/wl 
near the ground and entrainment zone follows a natural progression (increase, peak, and 
decrease) when transitioning from small-scale (2 km, A/^j ~ 1) to large-scale (30 km, A/^j ~ 
18) heterogeneity (Figure 12). Lateral variances {{v"^)/w1) are quite similar for all cases 
(not shown). 
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The increase in horizontal velocity fluctuations, at both the ground and Zi, is consis- 
tent with the presence of 2-D organized circulations induced by the strip-like heterogeneity. 
Hadfield et al. (1991) and Avissar and Schmidt (1998) both showed velocity fluctuations to 
increase at these locations in the presence of heterogeneity. However, Hadfield et al.’s (1991) 
largest heterogeneity scale (X/zi) was about 4 and they were therefore unable to trace the 
pattern back towards homogeneity with increasing X/zi as depicted in Figure 12. Avissar 
and Schmidt’s (1998) results, however, suggested that for this surface heat flux across the 
patches (~ 0.1 m K s“^) the fluctuations should trace the pattern like that in Figure 12. 

With the exception of the increase of {u"^)/wl near Zi, the response of the velocity 
variances to heterogeneity (Figures 12 and 13) is strikingly similar to the expected change 
between shear- and buoyancy-driven PBLs. Moeng and Sullivan (1994) showed (u"^)/wl 
to increase appreciably near the ground in the presence of a mean velocity shear. They 
also showed {w"^)fw'^ to decrease in the mid-PBL for increasing shear. Therefore at least 
with respect to velocity variances, these free convective PBLs in the presence of land-surface 
heterogeneity are somewhat analogous to free convective PBLs over homogeneous surfaces 
but with a mixed shear-buoyant forcing. 

6.3.2 Vertical Scalar Flux 

[Figure 14 about here.] 

For the DP PBLs, profiles of normalized vertical mixing ratio flux show only a small 
influence of surface heterogeneity (left panel. Figure 14). WP cases also show little influence 
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of heterogeneity (right panel, Figure 14), but the horizontally homogeneous cases reveal the 
influence of the soil partitioning of the available solar energy. As discussed in Section 6.2.1, 
case WP-DS partitions less available energy into latent heat and therefore more energy goes 
into sensible heat. Case WP-DS also entrains nearly three times more dry air than does 
WP-AS due to the much larger moisture jump across Zj. 

[Figure 15 about here.] 

Normalized profiles of vertical virtual potential temperature flux are linear through the 
PEL suggesting that our averaging is sufficient to assume quasi-steady conditions (Figure 
15). These results are counter to earlier findings (Avissar and Schmidt, 1998) who suggested 
that vertical potential temperature flux profiles are dramatically modified by heterogeneity. 
The profiles reveal very little modification due to land-surface heterogeneity no matter if the 
PEL is initially dry or wet (Figure 15). These results are consistent with the entrainment rate 
calculations presented in Figure 11, and suggest that the moisture component is truly not 
passive. In all cases, the profiles have a minimum flux (just below Z{) to surface flux ratio of 
about -0.12 which is smaller than the traditional value of -0.2 (Deardorff, 1979). According to 
Lilly (2002), this smaller minimum flux to surface flux ratio results from horizontal averaging 
across an undulating interface. 

7. Summary and Conclusions 

We have developed a large-eddy simulation code coupled with a fully interactive land-surface 
model. Initial temperature and moisture conditions are specified for the PEL and the soil 
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and are taken from SGP97 measurements. No mean wind is imposed such that the incoming 
solar radiation is the only external forcing in the simulations. The results show that soil 
moisture content is a crucial variable influencing PEL motions. 

Visualizations of instantaneous flelds reveal the signature of cellular convection at the 
surface and in the soil. Instantaneous slices of the PEL show dramatic modification of the 
PEL structure resulting from soil moisture heterogeneity. Instantaneous surface fluxes also 
depict the dynamic response of the forcing to the PEL. 

Consistent with previous research, under free convective conditions, heterogeneity gener- 
ates organized motions. Heterogeneity at scales ranging from \jzi ~ 4-9 induces the greatest 
PEL response. The organized motions induced at these scales of heterogeneity decrease the 
exchange coefficients for (momentum, scalars) by about (6, 9)% compared to the homoge- 
neous average soil cases (AS), and increase the bulk turbulence kinetic energy of the PEL 
by about 20%. The generation of these organized motions and the subsequent increase in 
turbulence kinetic energy, however, do not dramatically influence entrainment rates. 

Implications from these findings are that large-scale models that cannot resolve land- 
surface heterogeneity will slightly underestimate the surface fluxes of heat and moisture, 
but should not mis-estimate PEL entrainment. With respect to surface-layer turbulence 
statistics, heterogeneity makes the buoyantly forced PELs appear more like mixed shear- 
buoyancy PELs due to the induced organized motions. 
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Appendix A: Numerics 

For all variables the LES uses pseudo-spectral differencing in the horizontal. For velocities, 
the LES uses second-order centered-in-space differencing, while for scalars the monotone 
advection scheme of Koren (1993) is used. Spalart et al.’s (1991) third-order Runge-Kutta 
scheme advances all fields in time using a fixed Courant- Friedrichs- Lewy (CFL) number of 
0.63 for a complete 3-stage time step. Periodic boundary conditions are imposed in the hor- 
izontal. To allow the propagation of gravity waves out of the top of the domain, the upper 
boundary condition is that of Klemp and Durran (1983). Because the LES uses spectral dif- 
ferencing in the horizontal, the MPI implementation breaks up the domain in the vertical into 
slabs of horizontal data residing on a single processor. To perform the spectral decomposition, 
we use NCAR’s FFTPACK version 4.1 (http://www.scd.ucar.edu/softlib/FFTPACK.html). 

In the LSM, upwind finite differences are used for vertical derivatives and the Crank- 
Nicholson scheme is used for time differencing. At the lower boundary, the hydraulic diffu- 
sivity is assumed to be zero so that the soil water flux is is determined solely by gravitational 
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conduction. The soil heat flux is computed from a vertical temperature gradient between 
the lowest LSM model level and a specified soil temperature at an imaginary grid point one 
meter below the bottom of the model. Further details can be found in Mahrt and Pan (1984), 
Pan and Mahrt (1987), and Chang et al. (1999). To improve scalability when solving the 
LSM for the coupled surface boundary condition, we use MPI to break up the calculation in 
the x-direction. 
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Figure 1: Time evolution of the four vertical levels of soil moisture and temperature from the 1-D off-line 
LSM run initialized from measurements from Little Washita at 11:30 LST on Day 193 of SGP97. Initial soil 
conditions for the LSM for the wet and dry coupled simulations are taken at noon seven days apart. 
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Normalized Velocity Variances 


Figure 3: Comparison of velocity variances (u = (w"^ + |e}/tyj, w = + ~e)lwl, and SGS = 

derived from the current coupled LES-LSM code (solid lines) and those measured in the tank experiments 
of Willis and Deardorff (1974) case SI (x,#) and the LES results of Schmidt and Schumann (1989) (dashed 
line). Note that the bold typeface does not imply a vector quantity. The SGS energy (e) predicted by the 
LES is assumed isotropic such that one-third is added to each of velocity variances. 
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Figure 4: Comparison of horizontal velocity two-point correlations versus increasing lag in the 

lateral direction (A^) for three cases that are identical to case DP-15S except for the lateral extent of the 
domain. The second panel depicts vertical profiles of the square root of the horizontally- and time-averaged 
horizontal velocity variance normalized by for the same three cases. 
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Figure 5: Instantaneous (x,j/) slices of surface fluxes (sensible, latent and soil heat fluxes) and the corre- 
sponding soil temperature and moisture at the uppermost vertical level in the LSM (z = - 2.5 cm, which 
represents the vertical average over the top five centimeters in the soil) from DP- AS at the last time step in 
the simulation. 
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Figure 6: Instantaneous two dimensional slices of vertical velocity and water vapor mixing ratio at the first 
time step used in the averaging for case DP- AS (top two panels). The third panel depicts the coincident 
surface forcing (H, LE, -G) determined via the LES-LSM coupling. The lowest panel is the volumetric soil 
moisture at the uppermost soil level, z = - 2.5 cm. 
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X / Zi ^ / Zi 

Figure 9: Boundary-layer and time averaged turbulence kinetic energy normalized by versus the normal- 
ized heterogeneity wavelength (X/zi) for both the DP and WP cases. The line with symbols represents the 
heterogeneous soil cases. The solid, dashed and dotted lines with no symbols are the homogeneous soil cases. 
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A / zj A / zj 

Figure 11: Normalized entrainment rates {we/w^) versus the normalized heterogeneity wavelength {\/zi) 
for both the DP and WP cases. The line with symbols represents the heterogeneous soil cases. The solid, 
dashed and dotted lines with no symbols are the homogeneous soil cases. 
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Figure 12: Vertical profiles of normalized horizontal velocity variance, ( 1 (e) ) , for both the DP and 

WP cases. Lines by themselves represent heterogeneous cases. Lines with symbols represent homogeneous 
cases. The variable e is the subgrid-scale energy as predicted by the LES. 
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Figure 13: Vertical profiles of normalized vertical velocity variance, {{w^^) + for both the DP and 

WP cases. Lines by themselves represent heterogeneous cases. Lines with symbols represent homogeneous 
cases. 
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Figure 14: Vertical profiles of normalized vertical flux of water vapor mixing ratio ~\-r-u;q)/{w^Q^)) for 

both the DP and WP cases. Lines by themselves represent heterogeneous cases. Lines with symbols represent 
homogeneous cases. Tyjq is the subgrid-scale vertical mixing ratio flux predicted by the subgrid-scale model. 
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Figure 15: Vertical profiles of normalized vertical flux of virtual potential temperature {(w'^Oy -f 

for both the DP and WP cases. Lines by themselves represent heterogeneous cases. Lines 
with symbols represent homogeneous cases, is the subgrid-scale vertical virtual potential temperature 
flux predicted by the subgrid-scale model. 


49 





List of Tables 


1 Bulk properties for all cases. The quantities are: Zi, the PBL height; X/zi, 
the ratio of the scale of the heterogeneity (A) to the PBL height; w^., the 
convective velocity; H, the surface sensible heat flux (= pCpW^O^)] LE, the 
surface latent heat flux (= pLw^Q^)] G, the soil heat flux; the total 
surface buoyancy flux (= (w"0")+O.61 $o {w"q")); P, the Bowen ratio (H/LE). 

These bulk quantities are averaged horizontally over all space and over the 
entire averaging period. The angle brackets denoting this averaging process 
are omitted for convenience, n/a stands for not applicable since A is undefined 
for the homogeneous cases 51 


50 







» 


Case 

(km) 

A/^i 

(m s“^) 

H 

(W m-2) 

to 

G 

(W m-2) 

By, 
(m K 

P 

DP-WS 

1.66 

n/a 

2.35 

259.42 

138.97 

80.98 

0.23 

1.87 

DP-DS 

1.85 

n/a 

2.69 

358.74 

25.32 

78.22 

0.32 

14.17 

DP-AS 

1.79 

n/a 

2.58 

326.70 

59.91 

81.26 

0.29 

5.45 

DP-30S 

1.71 

17.55 

2.52 

311.00 

83.41 

77.82 

0.27 

3.73 

DP-15S 

1.75 

8.56 

2.54 

312.61 

84.30 

76.46 

0.28 

3.71 

DP-7.5S 

1.73 

4.33 

2.52 

310.51 

83.94 

78.59 

0.27 

3.70 

DP-5S 

1.72 

2.91 

2.51 

308.05 

82.97 

80.90 

0.27 

3.71 

DP-3S 

1.74 

1.72 

2.51 

307.14 

82.56 

80.94 

0.27 

3.72 

DP-2S 

1.74 

1.15 

2.51 

308.05 

82.86 

79.83 

0.27 

3.72 

WP-WS 

1.67 

n/a 

2.37 

264.72 

133.73 

80.68 

0.23 

1.98 

WP-DS 

1.79 

n/a 

2.67 

361.68 

23.77 

77.20 

0.32 

15.22 

WP-AS 

1.76 

n/a 

2.58 

329.94 

56.80 

81.12 

0.29 

5.81 

WP-30S 

1.70 

17.61 

2.52 

315.22 

79.73 

76.92 

0.28 

3.95 

WP-15S 

1.70 

8.82 

2.52 

316.11 

80.53 

75.96 

0.28 

3.93 

WP-7.5S 

1.72 

4.37 

2.53 

315.62 

80.78 

76.69 

0.28 

3.91 

WP-5S 

1.72 

2.91 

2.52 

312.61 

79.76 

79.71 

0.28 

3.92 

WP-3S 

1.72 

1.74 

2.52 

311.43 

78.88 

80.48 

0.28 

3.95 

WP-2S 

1.73 

1.55 

2.52 

312.05 

79.02 

79.73 

0.28 

3.95 


Table 1: Bulk properties for all cases. The quantities are: Zi, the PBL height; A/zj, the ratio of the scale 
of the heterogeneity (A) to the PBL height; tt;*, the convective velocity; H, the surface sensible heat flux 
(= pCpW^O^); LE^ the surface latent heat flux (— pLw^O^); G, the soil heat flux; the total surface 
buoyancy flux (= {w"0^') +0.61^o (iy"g")); /?, the Bowen ratio (H/LE), These bulk quantities are averaged 
horizontally over all space and over the entire averaging period. The angle brackets denoting this averaging 
process are omitted for convenience, n/a stands for not applicable since A is undefined for the homogeneous 
cases. 
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